
library(tidyverse)
library(survey)
library(here)

#============================================================#
# MODELS
#============================================================#

sub_dat <- read_rds(here("data","dat.rds"))
des <- svydesign(ids = ~PSU, strata = ~Strata, weights = ~MergeWgt, data = sub_dat, nest = TRUE)

#--------------#
# table 1
#--------------#
summary(svyglm(out ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))
summary(svyglm(out ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des, family = quasibinomial()))
summary(svyglm(out_cont ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))

#--------------#
# table a.1
#--------------#
summary(svyglm(out ~ treated_full * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))
summary(svyglm(out ~ treated_full * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des, family = quasibinomial()))
summary(svyglm(out_cont ~ treated_full * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))

#--------------#
# table a.2
#--------------#
summary(svyglm(out_age ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))
summary(svyglm(out_age ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des, family = quasibinomial()))
summary(svyglm(out_age_cont ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))

#--------------#
# table a.3
#--------------#
summary(svyglm(out ~ treated * newspaper + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))
summary(svyglm(out ~ treated * newspaper + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des, family = quasibinomial()))
summary(svyglm(out_cont ~ treated * newspaper + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))

#--------------#
# table a.4
#--------------#
summary(svyglm(out ~ treated * television + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))
summary(svyglm(out ~ treated * television + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des, family = quasibinomial()))
summary(svyglm(out_cont ~ treated * television + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))

#--------------#
# table a.5
#--------------#
summary(svyglm(out ~ treated_placebo * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag180, design = des))

#--------------#
# table a.6
#--------------#
summary(svyglm(out ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag3, design = des))
summary(svyglm(out ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag3, design = des, family = quasibinomial()))
summary(svyglm(out_cont ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag3, design = des))

#--------------#
# table a.7
#--------------#
summary(svyglm(out ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag30, design = des))
summary(svyglm(out ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag30, design = des, family = quasibinomial()))
summary(svyglm(out_cont ~ treated * radio + Ethnicity + Education + Income + Age + I(Age^2) + Gender + nattacks_lag30, design = des))

#============================================================#
# FIGURES
#============================================================#

#--------------#
# figure 1
#--------------#
read_csv(here("data","prop_dat.csv")) %>%
  ggplot(aes(x = date, y = interview, fill=radio_char)) +
  geom_bar(position="stack", width = 1, stat="identity", color=NA) +
  guides(fill=guide_legend(title="Radio")) +
  xlab("Date") +
  ylab("Respondents") +
  annotate("text", x = as.Date("April 25 2018",format = "%B %d %Y"), y= 300, label = "Taliban Radio\nBroadcast", size = 5, colour = "white") +
  annotate("text", x = as.Date("March 6 2018", format = "%B %d %Y"), y= 635, label = "Pre-Treatment", size = 5, colour = "white") +
  annotate("text", x = as.Date("May 22 2018",  format = "%B %d %Y"), y= 635, label = "Post-Treatment", size = 5, colour = "white") +
  theme(plot.title = element_text(hjust = 0.5),
        panel.grid.major = element_blank(),
        legend.position = "top"
  ) +
  scale_y_continuous(expand = expansion(mult = c(0.00001, 0)), limits = c(0,660)) +
  annotate("text", x = as.Date("April 25 2018",format = "%B %d %Y"), y= 300, label = "Taliban Radio\nBroadcast", size = 5) +
  annotate("text", x = as.Date("March 6 2018", format = "%B %d %Y"), y= 635, label = "Pre-Treatment", size = 5) +
  annotate("text", x = as.Date("May 22 2018",  format = "%B %d %Y"), y= 635, label = "Post-Treatment", size = 5) +
  geom_vline(xintercept = as.Date("May 5 2018",format = "%B %d %Y"), linetype="dotted")

#--------------#
# figure a.2
#--------------#
read_csv(here("data","diff_dat.csv")) %>%
  ggplot(aes(x = Wave, y = mean_q177, color = `Radio?`)) +
  geom_line() +
  ggtitle("Proportion of Respondents that Say the Taliban is Getting Stronger") +
  xlab("Survey Wave") +
  ylab("Proportion") +
  geom_vline(xintercept=c(39.5), linetype="dotted") +
  theme(plot.title = element_text(hjust = 0.5)
        ,text=element_text(size=19))

#--------------#
# figure a.3
#--------------#
read_rds(here("data","afgGED.RDS")) %>%
  ggplot(aes(x = date, y = n_attacks)) + 
  geom_point(color='black', shape=1) +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  xlab("Date") +
  ylab("Number of Attacks") +
  ggtitle("Taliban Attacks in Afghanistan Surrounding Survey Waves") +
  theme(panel.grid.minor = element_blank(),
        axis.text    = element_text(size=16),
        axis.title   = element_text(size=20),
        legend.text  = element_text(size=16),
        legend.title = element_text(size=16),
        plot.title   = element_text(hjust = 0.5, size=22)) +
  scale_y_continuous(expand = c(0, 0)) +
  annotate(geom = "rect", 
           xmin = as.Date("February 25 2018",format = "%B %d %Y"),
           xmax = as.Date("March 10 2018",format = "%B %d %Y"), 
           ymin = 0,
           ymax = 25,
           fill = "orange", 
           colour = "black",
           alpha = 0.3) +
  annotate(geom = "rect", 
           xmin = as.Date("May 15 2018",format = "%B %d %Y"),
           xmax = as.Date("May 29 2018",format = "%B %d %Y"),
           ymin = 0,
           ymax = 25,
           fill = "orange", 
           colour = "black",
           alpha = 0.3) 

#--------------#
# figure a.4
#--------------#
read_rds(here("data","sigar.RDS")) %>%
  ggplot(aes(x = year, y = enemy_attacks, fill=Color)) +
  geom_bar(position="stack", width = .175, stat="identity") +
  theme(axis.text    = element_text(size=16),
        axis.title   = element_text(size=20),
        legend.text  = element_text(size=16),
        legend.title = element_text(size=16),
        plot.title   = element_text(size=22)) +
  guides(fill=guide_legend(title="")) +
  scale_fill_discrete(breaks = "Study Year") +
  xlab("Year") +
  ylab("Number of Attacks") +
  ggtitle("SIGAR Enemy-Initiated Attacks Surrounding Survey Waves") +
  theme(axis.text    = element_text(size=16),
        axis.title   = element_text(size=20),
        legend.text  = element_text(size=16),
        legend.title = element_text(size=16),
        plot.title   = element_text(size=22),
        legend.position = "top")

#============================================================#








